home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / interval.lha / interval / interval.c next >
C/C++ Source or Header  |  1989-10-31  |  5KB  |  193 lines
  1. // This software is in the public domain.
  2. // Permission is granted to use this software for any
  3. // purpose commercial or non-commercial.
  4. // The implementer/distributer assume no liability for
  5. // any damages, incidental, consequential or otherwise,
  6. // arising from the use of this software.
  7.  
  8. // Implementation of class interval
  9. //
  10. // Most operators are trivial and could be inlined.
  11. //
  12. // Multiply is more interesting.  The complicated testing
  13. // may be more time consuming than the floating-point
  14. // operations it attempts to avoid.
  15. // How could this testing be streamlined?
  16. //
  17. // Note that divide blurts out error/warning messages using fprintf :-(
  18. // Strictly speaking, division is undefined if the divisor interval
  19. // contains zero.  This code is sloppy in allowing it, but it warns.
  20. //
  21. // No attempt is made to control rounding so accuracy is dubious.
  22. //
  23. // No attempt is made to detect/recover from overflow or underflow.
  24. // This is a serious deficiency.  Is there a 'portable' around it?
  25.  
  26. #include <CC/stdio.h>
  27. #include "interval.h"
  28.  
  29. interval interval::operator+(interval I)
  30. {
  31.     interval temp(lo_bound, hi_bound);
  32.     return temp += I;
  33. }
  34.  
  35. interval interval::operator-(interval I)
  36. {
  37.     interval temp(lo_bound, hi_bound);
  38.     return temp -= I;
  39. }
  40.  
  41. interval interval::operator*(interval I)
  42. {
  43.     interval temp(lo_bound, hi_bound);
  44.     return temp *= I;
  45. }
  46.  
  47. interval interval::operator/(interval I)
  48. {
  49.     interval temp(lo_bound, hi_bound);
  50.     return temp /= I;
  51. }
  52.  
  53. interval &interval::operator+=(interval I)
  54. {
  55.     lo_bound += I.lo_bound;
  56.     hi_bound += I.hi_bound;
  57.     return *this;
  58. }
  59.  
  60. interval &interval::operator-=(interval I)
  61. {
  62.     lo_bound -= I.lo_bound;
  63.     hi_bound -= I.hi_bound;
  64.     return *this;
  65. }
  66.  
  67. // Multiply is a bit complicated.
  68. // A naive version simply takes the minimum of all combinations
  69. // as the new lo_bound and the maximum as the new hi_bound.
  70. // But this involves four double precision floating multiples.
  71. // A quick examination of the signs of the intervals reveals
  72. // that there are nine possible cases of which only one
  73. // requires all four multiplications.  The other cases only
  74. // require two multiplications.  We assume that it is faster to
  75. // compute and dispatch to the right case than to perform
  76. // the extra multiplies.  If not, use the naive scheme.
  77. //
  78. // The cases can be described as follows:
  79. // Each interval is either non-negative (lo bound >= 0),
  80. // non-positive (hi bound <= 0), or crosses(includes) zero,
  81. // including some numbers on either side.
  82. // If we label the intervals, A and B, corresponding to
  83. // "this" and "I", we get the following matrix:
  84. //
  85. //        Bcross        Bnonneg        Bnonpos
  86. //
  87. //    Across
  88. //
  89. //    Anonneg
  90. //
  91. //    Anonpos
  92.  
  93. interval &interval::operator*=(interval I)
  94. {
  95.     enum possibility {
  96.     AcrossBcross, AcrossBnn, AcrossBnp,
  97.     AnnBcross, AnnBnn, AnnBnp,
  98.     AnpBcross, AnpBnn, AnpBnp
  99.     } choice;
  100.  
  101.     if (lo_bound >= 0.0) choice = AnnBcross;
  102.     else if (hi_bound <= 0.0) choice = AnpBcross;
  103.     else choice = AcrossBcross;
  104.  
  105.     if (I.lo_bound >= 0.0) choice += 1;
  106.     else if (I.hi_bound <= 0.0) choice += 2;
  107.  
  108.     switch (choice)
  109.     {
  110.     case AcrossBcross: {
  111.             double  HL = hi_bound*I.lo_bound,
  112.                 HH = hi_bound*I.hi_bound,
  113.                 LL = lo_bound*I.lo_bound,
  114.                 LH = lo_bound*I.hi_bound;
  115.             lo_bound = HL<LH?HL:LH;
  116.             hi_bound = LL>HH?LL:HH;
  117.         }
  118.         break;
  119.  
  120.     case AcrossBnn:
  121.         lo_bound *= I.hi_bound;
  122.         hi_bound *= I.hi_bound;
  123.         break;
  124.  
  125.     case AcrossBnp: {
  126.         double new_hi_bound = lo_bound * I.lo_bound;
  127.  
  128.         lo_bound = hi_bound*I.lo_bound;
  129.         hi_bound = new_hi_bound;
  130.         }
  131.         break;
  132.  
  133.     case AnnBcross:
  134.         lo_bound = hi_bound*I.lo_bound;
  135.         hi_bound *= I.hi_bound;
  136.         break;
  137.  
  138.     case AnnBnn:
  139.         lo_bound *= I.lo_bound;
  140.         hi_bound *= I.hi_bound;
  141.         break;
  142.  
  143.     case AnnBnp: {
  144.         double new_hi_bound = lo_bound * I.hi_bound;
  145.  
  146.         lo_bound = hi_bound*I.lo_bound;
  147.         hi_bound = new_hi_bound;
  148.         }
  149.         break;
  150.  
  151.     case AnpBcross:
  152.         hi_bound = lo_bound*I.lo_bound;
  153.         lo_bound *= I.hi_bound;
  154.         break;
  155.  
  156.     case AnpBnn:
  157.         lo_bound *= I.hi_bound;
  158.         hi_bound *= I.lo_bound;
  159.         break;
  160.  
  161.     case AnpBnp: {
  162.         double new_hi_bound = lo_bound * I.lo_bound;
  163.  
  164.         lo_bound = hi_bound*I.hi_bound;
  165.         hi_bound = new_hi_bound;
  166.         }
  167.         break;
  168.  
  169.     default:
  170.         fprintf(stderr, 
  171.             "Bad interval (low bound > hi bound) [%f, %f] * [%f, %f]\n",
  172.             lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  173.         lo_bound = hi_bound = 0.0;
  174.         break;
  175.     }
  176.  
  177.     return *this;
  178. }
  179.  
  180. interval &interval::operator/=(interval I)
  181. {
  182.     if (I.lo_bound == 0.0 && I.hi_bound == 0)
  183.     fprintf(stderr, "Error: Division by zero attempted - [%f,%f]/[f,%f]\n",
  184.         lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  185.     else {
  186.     if (I.lo_bound < 0 && I.hi_bound > 0)
  187.         fprintf(stderr, "Warning: Division by zero possible - [%f,%f]/[%f,%f]\n",
  188.         lo_bound, hi_bound, I.lo_bound, I.hi_bound);
  189.     interval inv(1.0/I.lo_bound, 1.0/I.hi_bound);
  190.     return inv *= *this;
  191.     }
  192. }
  193.